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We study the ground-state Wigner bilayers of pointlike particles with Yukawa pairwise inter¬ 
actions, confined to the surface of two parallel hard walls at dimensionless distance r]. The model 
involves as limiting cases the unscreened Coulomb potential and hard spheres. The phase diagram of 
Yukawa particles, studied numerically by Messina and Lowen [Phys. Rev. Lett. 91 (2003) 146101], 
exhibits five different staggered phases as r) varies from 0 to intermediate values. We present a 
lattice summation method using the generalized Misra functions which permits us to calculate the 
energy per particle of the phases with a precision much higher than usual in computer simulations. 
This allows us to address some tiny details of the phase diagram. Going from the hexagonal phase 
I to phase II is shown to occur at ?7 = 0, which resolves a longtime controversy. We find a tricritical 
point where Messina and Lowen suggested a coexistence domain of several phases which was sug¬ 
gested to divide the staggered rhombic phase into two separate regions. Our calculations reveal one 
continuous region for this rhombic phase with a very narrow connecting channel. Further we show 
that all second-order phase transitions are of mean-field type. We also derive the asymptotic shape 
of critical lines close to the Coulomb and hard-spheres limits. In and close to the hard-spheres limit, 
the dependence of the internal parameters of the present phases on rj is determined exactly. 

PACS numbers: 68.65.Ac, 52.27.Lw, 82.70.Dd 


I. INTRODUCTION 


Most organic or anorganic surfaces of mesoscopic ob¬ 
jects (macromolecules or colloids) become charged when 
immersed in polar solvents such as water. These sol¬ 
vents provide favorable environments for free charges 
(“counter-ions” to the charged surfaces) which interme¬ 
diate an effective interaction among the mesoscopic ob¬ 
jects. At low temperatures, and in particular at T = 0 
when the system is in its ground state, counter-ions be¬ 
tween two charged plates crystallize into bilayer Wigner 
structures which are important in understanding anoma¬ 
lous phenomena such as like-charge attraction or over¬ 
charging Bilayer Wigner crystals describe several 

real physical systems in condensed and soft matter, such 
as semiconductors Q, quantum dots Q and dusty plas¬ 
mas 0. Confined systems of charged colloidal particles 
were reviewed recently Q, both from experimental and 
theoretical point of view. 

From the particle models studied in this paper, we 
start with the neutral Coulomb system of say elementary 
pointlike charges —e with 1/r interaction between two 
parallel plates of the same homogeneous surface charge 
density ae at distance d, the phase diagram at T = 0 
depends on a single dimensionless parameter rj = dy/a. 
According to the Earnshaw theorem 0, particles will 
stick symmetrically on the surface of the plates. Five 
distinct phases were detected to be stable, i.e. provid¬ 
ing global minimum of the energy, as ry is changing from 
0 to oo [IMl. The lattice structures are the same on 
both plates and they are shifted laterally with respect 
to one another. Structures I, III and V are rigid (Fig. 
[1]), i.e. they have fixed (ry-independent) primary cells. 
Structures II and IV are soft (Fig. [2]), the shape of their 
primary cells is varying with ry. 
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FIG. 1. Rigid structures I, III and V of particles on two 
parallel plates; open and filled symbols correspond to particle 
positions on the opposite layers. 


Structures I, II and III correspond to the staggered 
rectangular lattice, see Fig. [^left. The primitive trans¬ 
lation vectors of Bravais lattice are 


ai = a(l, 0), 02 = a(0. A), 


1 



( 1 ) 


The lattice spacing a within one layer is determined by 
the electroneutrality condition. Two rectangular struc- 
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FIG. 2. Soft structures II and IV. 


tures, one within each layer, are shifted with respect to 
each other by the vector 

c = 0 ( 01 + 02 ) (2) 

with a = 1/2. The rigid structure I has the aspect ratio 
A = and it arises for ry = 0 because the two layers 
merge into a Wigner monolayer which is known to be 
hexagonal (or, equivalently, equilateral triangular) la- 
structure III consists of a square lattice with A = 1. 
Phase II with 1/3 + 1 interpolates continuously 

between structures I and III. 

Phase IV consists of two staggered rhombus lattices 
(Fig. [2] right). One rhombus structure has the angle (j) 
between the primitive translation vectors 

Oi=a(l,0), 02 = a(cos(/, sini/), a = —= , (3) 

V a sin 0 

In general, this phase has two variants according to the 
lateral shift (I5|) between the opposite sublattices. The 
version IVA has a = 1/2 whereas for IVB the shift pa¬ 
rameter 1/3 < a < 1/2. For Coulomb bilayers, only 
phase IVA takes place. 

Phase V corresponds to two shifted hexagonal lattices. 
In a single layer, the elementary cell is the rhombus with 
the primitive translation vectors 

ai=a(l,0), a2 = ^(l,\/3), a = (4) 

The lateral shift between the opposite lattices c is given 
by (m with a = 1/3. 

The transitions between phases II —>■ III and III —>■ IV 
are continuous (of second order), while the transition 
IV —?> V is discontinuous (of first order). In order to de¬ 
scribe these phase transitions, a new analytic approach 
to Coulomb bilayers was proposed in Ref. [l^. Using 
a series of transformations with Jacobi theta functions, 
the energy of the five phases was expressed as series of 
generalized Misra functions which converge very quickly. 
Near critical points, the generalized Misra functions can 
be expanded easily in powers of the order parameter and 
the corresponding energies posses the Ginsburg-Landau 


form. This allows one to specify the critical points with 
an arbitrary prescribed accuracy and to derive the mean- 
field critical behavior of the order parameter. Also the 
existence of phase I at ry = 0 only was confirmed. This re¬ 
sult was in contradiction with numerical approaches like 
Ewald technique [l^ and Monte Carlo simulations [l^ 
which predicted an extremely small, but finite, stability 
interval of ?y’s for phase I. 

Colloidal particles or particles in highly charged dusty 
plasmas usually interact via Yukawa potential (l8l | due to 
the Coulomb potential screening by additional microions 
in the system. The Yukawa pair potential of particles at 
distance r is defined by 


-ftr 

U(r) = Uo -, (5) 

KT 

where k is the inverse screening length and the ampli¬ 
tude Vq = Z^Kexp (Ki?)/e(l + kR)'^, with Z being the 
charge of one particle and e ~ eo is the permittivity for 
dusty plasma. When k is large, R becomes the radius 
of a hard sphere as V{r) oc exp [^(i? — r)] is exponen¬ 
tially large for r < R and negligible otherwise. The re¬ 
lation Vq oc K keeps the limit k —>■ 0 of Eq. m finite, 
yielding the proper Coulomb formula. Thus the limiting 
cases K —>■ 0 and k ^ 00 correspond to the unscreened 
Coulomb and hard-spheres interaction potentials, respec¬ 
tively. We shall work in units of Vq = 1 ■ For two parallel 
plates at distance d, the phase diagram depends on two 
dimensionless parameters 

ry = ^/ad, A = Kd. (6) 


A system of hard-sphere particles between two parallel 
hard plates was studied by computer simulations in the 
past II 9 I - I 22 I : numerical methods were reviewed recently 
in Ref. [^. For small values of ?y, the ground-state 
crystal structures involve Wigner bilayers I-V, includ¬ 
ing phase IVB with two varying parameters 0 and a. 
For large values of ry, phase-V bilayer transforms itself to 
crystalline multilayers, with particles entering the region 
between the plates, such as multiple square and hexago¬ 


nal layers [19|, rhombic [20[ and prism superlattices |21| . 

A similar phase diagram was obtained for the general 
Yukawa potential. For small values of ry, although Earn- 
shaw theorem does not apply to Yukawa particles, 
the particles stick symmetrically to plates and with in- 
creasing r? they constitute successively Wigner I-V bi¬ 
layers 


24| . In the region of large values of ry, in close 
analogy with confined hard spheres, some of the parti¬ 
cles will move in the interior of the domain between the 
plates and create multilayers [ 2 ^. 

In this paper, we shall concentrate on Wigner bilay¬ 
ers of pointlike particles interacting via Yukawa poten¬ 
tial. The original numerical work of Messina and Lowen 


24j determined the phase diagram of the Yukawa sys¬ 


tem which exhibits single and double reentrant transi¬ 
tion. We shall apply a straightforward extension of the 
recent analytic method [13 which provides us with high 
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precision calculations to shed more light on important 
tiny details of the phase diagram. For any A, the tran¬ 
sition from phase I to II is shown to occur directly at 
77 = 0, which solves a longtime controversy. We re¬ 
call that this scenario was anticipated only in the hard- 
spheres limit A —>■ 00 M- We find a tricritical point 
where Messina and Lowen suggested a coexistence do¬ 
main of several phases which should divide one staggered 
rhombic phase into two separate regions. Our calcula¬ 
tions reveal one continuous region for this rhombic phase 
with a narrow connecting channel. Closed-form formu¬ 
las for critical lines between various phases, expressed in 
terms of generalized Misra functions, permit us to deter¬ 
mine the asymptotic Coulomb A —>■ 0 and hard spheres 
A —>■ 00 shapes of these lines. The expansions of the struc¬ 
ture energies around second-order transition points and 
the determination of the order parameter can be done 
analytically, which enables us to derive the critical be¬ 
havior of the Ginzburg-Landau type. In and close to the 
hard-spheres limit, the 77 -dependence of the internal pa¬ 
rameters of the phases is determined exactly. 

The paper is organized as follows. In Sec. im we de¬ 
rive the expression for the energy per particle of phase II 
(phases I and III being its special cases) in terms of the 
generalized Misra functions. The fact that going from 
phase I to II occurs at 77 = 0 is shown in Sec. IIIII The 
second-order transition between phase II and III and the 
corresponding mean-field critical behavior are described 
in detail in Sec. El The expression for the energy of 
phase IVB (with phases IVA and V as its special cases) 
is derived in Sec. lYl The second-order transition be¬ 
tween phases III and IVA is described in Sec. |E The 
first-order transitions between phases IVA-V, IVA-IVB 
and IVB-V are discussed in Sec. I VIII The dependence 
of the energy on the dimensionless distance 77 , for fixed 
values of A, is the subject of Sec. IVIIIl The ? 7 -dependence 
of the internal structure parameters of phases present in 
and near the hard-spheres limit is derived in Sec. ESec. 
1x1 is the Conclusion. Auxiliary formulas for the general¬ 
ized Misra functions and for the critical lines are given 
in Appendices A-D. 


II. ENERGY OF STRUCTURES I, II AND III 

We aim at deriving the interaction energy per par¬ 
ticle Ell for the structure II with the aspect ratio A, 
phases I and III being its special cases with A = 
and A = 1, respectively. The energy consists of two 
parts: the intralayer energy i?intra sums the contribu¬ 
tions from all particles in the same layer as the reference 
one while the interlayer energy Ai„ter involves all parti¬ 
cles from the opposite layer. To express the energy per 
particle as a quickly convergent series, we shall apply a 
three-step method from Ref. m. 

The occupied lattice sites within one layer are num¬ 
bered as r = jai + ka 2 where the primitive vectors ai 
and 02 are dehned in m and j, k run over all integers. 


except for the reference site (0,0). The intralayer inter¬ 
action of a reference particle is thus given by 


Ai 


intra 


=) i: 


(j,/c)#(0.0) 


exp 

Ka\J p + fc^A^ 


(7) 


To evaluate lattice sums of Yukawa potentials, we shall 
often use the integral representation (see e.g. 


1 


KT 




( 8 ) 


The intralayer energy per particle is then expressible as 


T/intra — 


laKy/i: Jq y/l 


^ dt 

e “i* 


- 1 


dt 


2y/nX Jq y/t 


e 


E 

, j,k 


O 3 (e-*^)03 (e-^) -1 


(9) 


where we substituted tA —>• t and introduced the Jacobi 
theta function 03 (q, 0 ) = 03 (q) = 

The Wigner lattice on the opposite layer at distance 
d is shifted by the vector (oi -|- a2)f2. The square of 
the distance between the reference particle and the par¬ 
ticles on the opposite layer becomes = (j — l/2)^a^ -I- 
{k — 1/2)‘^PA^ + d^. Proceeding analogously as in the 
previous case, we get for the interlayer energy 


Aintpr — 


dt -X, — 


2-\/^A Jg y/t 


0 Arj^t 


O 2 (e-‘^) 02 


( 10 ) 


where another Jacobi theta function 02{q) = 
was introduced. 

The total energy per particle Eu is a sum Ajntra + 
Ainter- Using the Poisson summation formula 


^ Q-U+^ft = Y 7 E (11) 


it can be easily shown that in the limit t —>■ 0 the product 
of theta functions 0 m(e“*) 6 >m(e“‘) ~ 7 r/t for both m = 
2,3. In the unscreened Coulomb limit A —0, this would 
lead to the divergence of the corresponding integrals due 
to the lack of the neutralizing background charge. We 
“artificially” subtract the singular 7r/t terms from the 
products of theta functions and simultaneously add the 
same singular terms and integrate them explicitly, with 
the result 
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This corresponds to adding and subtracting the back¬ 
ground interaction energy 


TT At' / TT 


.0 477 -^ 


TT ^ _ j t ^^ 


e-/c^dA _ 1 ^ 

TT 




(13) 


dt_ ^ 


Vi 


0 477*2 


(«-“) »3 (. 




- 1 - 




The procedure is inevitable in the Coulomb A —>■ 0 limit. 
For a positive A > 0, the procedure is not necessary but 
it enhances substantially the convergence properties of 
the obtained series. 

The integration region [0, 00 ] in ([9]) can be split into 
intervals [0,7r] and [ 77 , 00 ]. Using the Poisson summation 
formula (HU, the integral over [tt, 00 ] can be rewritten as 



Similarly, 

dt 


L Vi 


e 




dt —^ 
_^e ‘‘1 

Vi 


O 2 (e-*^) 02 (e-^) - ■ 

04 (e“‘^) 04 (e“^) - 1 3 (15) 


where we introduced the Jacobi theta function 04 (q) = 

E,(-i)'/- 

Finally, in close analogy with Ref. M we apply once 
more the Poisson summation formula m for each term 
in the integration from [0, tt]. The final formula for the 
energy reads 


( 0 , A^/(47r^?7^) +fA) + Z 3/2 (O, A^/(47r^?7^) -h j^/A)] - ttzi/s ^0, ^^ 2 ^ 2 ^ 

00 

-f 2 ^(-l)^ [2:3/2 ( 7 ^^^ 7 ^AV( 47 ^^b^) + j^A) + Z3/2 {VV+ f /^)] - T^Zi/a (AV( 4 ??^), 0 ) 

i=i 

00 00 

+4 {-iy{-lVz 3 / 2 {VVV'^/{‘iT^‘^V)+f/A + VA)+4: ^ Z 3/2 {0,X'^/{4:VV) + f/A +VA) 

j,k—l j,k—l 

00 

+2^(-l)^' [2:3/2 {VV,X^/VT^^V)+f^)+Z3/2 (7^^^?^AV(47^^?7^) + jVA)] -7704/2 (AV( 47 ?^), ? 7 ^) 

7=1 

00 00 

4"^^ [^3/2 (V/{4V),VA) + Z 3/2 {V/ (4??^), j^/A)] J- 4 ^ Z 3/2 (A^/(4ry^), j^/A J- VA) 

j=l 7.fc=l 

00 2 

+4 ^ Z 3/2 [X^I{diV),V + (j - l/2)^/A + {k- l/2fA] I -h 77 (1 -h e"^) . (16) 

7A=1 ^ 


Here, we introduced the function 

= l -e-^V-y/y (17) 

It is a generalization of the well-known Misra function 
(^ . corresponding to x = 0, commonly used in lattice 
summations. The functions Zi,{x, y) with half-integer val¬ 
ues of v can be expressed in terms of the complementary 
error function, see Appendix A. This permits us to use 
very effectively the MATHEMATICA software and to de¬ 
rive in Appendix A their asymptotic forms for (x finite, 
y —>■ 00 ) and [y finite, x —>■ 00 ). The series in the gen¬ 
eralized Misra function dm) is quickly converging; for 
the known A = 0 Coulomb cases m , the truncation of 
the series over j,k at M = 1, 2, 3,4 reproduces the exact 


value of the energy up to 2, 5,10,17 decimal digits, re¬ 
spectively. This accuracy even improves itself for A > 0, 
so in our numerical calculations we keep the truncation 
of the series at M = 5. 

The formula (HH) is symmetric with respect to the 
transformation A —>■ 1/A. This symmetry corresponds 
to an obvious invariance of the energy with respect to the 
lattice rotation around one point by 90 degrees. 

III. GOING FROM PHASE I TO II 

As was mentioned in Introduction, numerical ap¬ 
proaches [H, M predicted that phase I has a region 
of stability [0, fy] with a very small rj > 0 and there is 
a second-order transition between phases I and II. This 
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small region was expected to vanish (iy = 0) in the hard- 
spheres limit A —>■ oo [l^. But in the paper [a it was 
shown both analytically and numerically that 77 = 0 in 
the unscreened Coulomb limit A —>■ 0 , i.e. phase I exists 
only for rj = 0 . There is no singularity in the ground- 
state energy, so going from phase I to phase II is not a 
phase transition in the usual sense. In what follows, we 
derive the same results for any positive A. 

We know that A = -x/S for phase I at ry = 0 . Let us 
assume that for 7y > 0 we have A = -^/S — e with a small 
e and, in close analogy with Ref. , expand the energy 
m in Taylor series: 

L;ii(V 3 - e, 7 y. A) = £^11(^3, ry. A) -b /i(?y, A)e 

+/2(»y, A)e^-b (!I(e^), ( 18 ) 

where the expansion functions /i(?y. A) and /2(ry, A) are 
written explicitly in terms of the generalized Misra func¬ 
tions in Appendix B. For given rj and A, the extremum 
of the energy (fT51) occurs at e* given by 


^£;ii(\/ 3 - e,?y. A) 
implying 


/i(7y. A) -b 2/2(77, A)e* = 0 , 
( 19 ) 


( 20 ) 


For the unscreened Coulomb case A = 0 it has been 
shown in that 

- A*(7y, 0 ) = = 7-14064 ... 77^ + 0(77^). 

This extremum is the minimum of i?ii(e). 

In the case of A > 0 , it is shown in Appendix B that 
for 77 <C A the coefficient functions can be approximated 
by 

/i( 77,A) «-^^^e ^-bO^^Ty^e 7^7%^, 

The extremum 

73 - A777, A) = = 27 + 0 ( 7 ) ( 23 ) 


interestingly does not depend in the leading order on 
A. It corresponds to the minimum of energy £'n(e) as 
d‘^Eii {'/3 — e, 77, A)|e=e. = 2/2(77, A) > 0 . For A = 1 and 
77 = 0 . 01 , we checked the result ( 1 ^ numerically in Fig. 
[21 One can see that Aii(e), calculated using the com¬ 
plete formula ( 1161 ) truncated at M = 5 has a minimum 
rather close to the value e* = 0.0002 predicted by our 
asymptotic formula (Uni- 

We conclude that phase I is stable only at 77 = 0 and for 
an arbitrarily small positive 77 we enter the region of phase 



FIG. 3. Eii[A,rf) — Ei{r]) as a function of — A for the 
fixed values of A = 1 and 77 = 0.01. The value of e* = 0.0002, 
which provides the energy minimum according to the asymp¬ 
totic formula (1231) is depicted by the vertical dashed line for 
comparison. Note that the energy differences are extremely 
small. 


II. It is interesting that the asymptotic 77 —>■ 0 predictions 
for the unscreened Coulomb A = 0 case m and for A > 0 
(1^21) exhibit the same rf' dependence, but there is a skip 
in the prefactors from 7 . 14064 ... at A = 0 to 2 for A > 0 . 
The fact that in the previous works [H, [H, [IJ phase I 
was detected also for small positive values of 77 is probably 
related to extremely small deviation of 73 — A* oc 77^ 
which are “invisible” by standard numerical methods. 


IV. SECOND-ORDER TRANSITION BETWEEN 
PHASES II AND III 


Let us parametrize A = exp(e). The symmetry A —>■ 
1 /A of the energy (fT21) is then equivalent to the trans¬ 
formation e —>■ — e and the energy is an even function of 
e. The Ginsburg-Landau form of its expansion around 
e = 0 reads as 


rj, A) = £1111(77, A) -b g 2 {v, A)e^ + Qiiv, A)e^ -b ... 

( 24 ) 

The explicit expression for <72 is given in Appendix C 
and a rather cumbersome expression for 34 is also at our 
disposal. The critical point is given by the vanishing of 
the prefactor 


32 ( 77 '=) = 0 . ( 25 ) 

We used this equation to get the (dashed) critical line 
between phases II and III in Fig. U) Our definition of 
77 differs from that of the dimensionless distance in the 
paper of Messina and Lowen namely 77^ = 77 ml- 
To maintain the full comparability, we shall present the 
phase diagram using the variable 77^. 
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FIG. 4. Phase diagram of the Yukawa bilayer. Dashed lines 
denote the second-order phase transitions, solid lines corre¬ 
spond to the first-order phase transitions. The important 
data for the hard-sphere limit A oo are added on the top. 



FIG. 5. Order parameter close to the critical point of the 
transition II-III for three values of A = 1,10,100. Full lines 
follow from numerical minimization of the energy (I16II . The 
slope of lines is close to /3mf = 1/2. Dashed lines represent 
the asymptotic t] relation (I28II . 


A. Critical behavior 


To obtain the critical behavior, we note that the func¬ 
tions g2 and 54 in Eq. ( 1241 ) behave in the vicinity of the 
critical point {rj^, A'^) as follows 

92iv, A) = 52i(A=)(r7^ - v) + 0[{r]^ - r])% 

54(?7, A) = 54o(A‘')- f - 77), ( 26 ) 

where g2i(A'^) < 0 and g4o{X‘^) > 0 for all A°. The mini¬ 
mum energy is reached at e* ~ A* — 1 given by 

^£;ii(e",?7, ~ 252(7?, A)e* -h 454(5, A)(e*)^ = 0 . 

( 27 ) 

For 5 > 5'^, there is only one solution e* = 0 which 
corresponds to the square lattice of phase III. For 5 < 5°, 
we get one trivial (unphysical) solution e* = 0 and two 
non-trivial conjugate solutions ±e* with 


e 


* 


32(77, A) 
234(7?, A) y 


32i(A^) 

2340(A°) 



- 3, 


( 28 ) 

5 —> ( 5 *^) . The order parameter e* oc is thus 

associated with the mean-field critical index /3 mf = 1/2 
for every A > 0. The dependence of A — 1 on 5 ° — 5 is 
shown in Fig. [S]for three values of A = 1,10,100. Near 
the critical point ( 5 “ — 5 small), the asymptotic relation 
( 1 ^ (dashed lines) fits perfectly the numerical data from 
minimization of the energy Eu (jl2p (full lines). In the 
logarithmic plot, for all values of A the slope of A — 1 
vs. 5 “^ — 5 is very close to 0.5 in the region of small and 
intermediate values of 5 '^ — 5 , confirming the value 1 /2 of 
the mean-field critical index /3mf for all values of A. 

From Eq. (IMl) . the energy difference of phases II and 


III close to the critical point is given by 

F„(e% 5, A) - Euiiv, A) ^ -^^( 3 ^ - 3 )". ( 29 ) 

The critical singularity should be of type (5'^ — 5)^““ 
implying the mean-field critical index omf = 0 for any 

A. 

To obtain another two critical indices, we add to the 
energy dH]) the symmetry-breaking term —he, where a 
small positive external field /i —>■ O'*' is linearly coupled 
to the order parameter. The optimization condition for 
the energy with respect to e now takes the form 

252(5, A)e* + 454(5, A)(e *)3 - h = 0 . ( 30 ) 


At the critical point, since 52(3^^, A'^) = 0 and 54(5"^, A'^) = 
540 (A'^), we find from (l 30 t that 


e 


* 


h 


1/3 


.4340(A"). 


( 31 ) 


This critical singularity should be of type which 

leads to the mean-field critical index Amf = 3 for any 
A. Performing the derivative of Eq. ( 1501 ) with respect to 
h, we find for the field succeptibility close to the critical 
point: 


de* 


dh 


h —0 


1 1 

-4521 (A°) 5^-5’ 


3 ^ ivl ■ ( 32 ) 


The corresponding critical singularity (3 —3°)”"*' leads to 
the mean-field critical index 7 mf = 1 for arbitrary A. 

It is easy to verify that our mean-field critical indices 

omf = 0 , / 3 mf = 2^ 7mf = 1 , ^mf = 3 ( 33 ) 
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fulfill two standard scaling relations 

2-a = 2/3 + 7 = /3(d + l). 


small again compared to the leading ones. In the last 
sum the Z 7 / 2 (-, •) term has zero prefactor for j = k. We 

(34) 

are left with the three-terms expression 


Since there are no fluctuations in our system at zero tem¬ 
perature, the critical indices rj and v, related to the par¬ 
ticle correlation function, are not defined. 


B. Coulomb A —0 limit of the critical line 

We reproduce 7/^(0) = 0.2627602682 [l^ in the 
Coulomb A —^ 0 limit. It is shown in Appendix C that 
the asymptotic A —>■ 0 shape of the critical line between 
phases II and III is parabolic: 

A"«iC23[7?"-77"(0)], C23 ~ 24.173744. (35) 

This formula is compared to the critical line evaluated 
numerically by using the relation (1251) in Fig. [S) 



FIG. 6. The critical line between phases II and III near the un¬ 
screened Coulomb A —>■ 0 limit. Full line follows from the nu¬ 
merical evaluation by using the relation < 72 ( 7 , A) = 0. Dashed 
line corresponds to the asymptotic formula (1351) . 


C. Hard-spheres X ^ 00 limit of the critical line 


52 ( 7 /, A) Z 7/2 1 ) - ^ 5/2 1 ) 

2) ’ ^ » 1- (36) 

Applying the asymptotic formula (IA9|) . we rewrite the 
rhs of this expression as 



A 

V 


V _ AV+ + 1 

A 7/2 + 1 


(37) 

The critical condition <72 ( 77 , A) = 0 implies a transcenden¬ 
tal formula for 77 (A): 


1+ . , 




(38) 


The exponential term can equal to the rational one only 
if I — •\/772 + 1/2 is close to zero, i. e. 77 ^ —1/2 in the 
A —>■ 00 limit as expected. The next terms of the large-A 
expansion of 77 (A) can be derived straightforwardly, with 
the result 


InA ln2 ,^/^ln^AA 


(39) 


In general, the series contains the terms of the form 
(lnA)"*/A"' where m, n are integers such that 0 < m < n. 
The first correction of type (lnA)/A explains a slow con¬ 
vergence of the results as A —>■ 00 . 

The asymptotic formula (l3^ . taken for 77 ^, is plotted 
in Fig. [7] by the dash-dotted line. We see that it repro¬ 
duces adequately the numerical results for the critical 
line (solid line) in a large region of the phase diagram. It 
can be shown that the next term of the series dMl) reads 
as 31n^ A/( 2 ^/^A^); plotting the asymptotic formula (IMl) 
with this term included makes the difference with the 
numerical solid line invisible by eye. 


In the hard-spheres limit A —>■ 00 , the critical point for 
the II —>■ III transition is (t/*^)^ —>■ 1/2 [^. The conver¬ 
gence to this value is extraordinarily slow. Let us analyze 
this limit in the critical equation < 72 ( 77 , A) = 0. Applying 
the asymptotic formulas for the generalized Misra func¬ 
tions (|A8I) and (IA9p to < 72 ( 77 , A) given by the series (|C1I) . 
most summands become exponentially small compared 
to the few leading terms proportional to exp(—A/ 77 ). In 
particular, we can neglect completely the first four sums 
in Eq. m since all terms behave as exp(—cA^) and the 
sixth sum because we get at least exp(— v^A/t/) for the 
j = k = 1 term. Those leading terms appear in the fifth 
sum with j = 1 and the seventh sum with j = k = 1. 
The ones with e. g. / = 2, fc = 1 etc. are exponentially 


V. ENERGY FOR STRUCTURES IVA, IVB 
AND V 

It was already mentioned that structures IVA, V and 
even III are special cases of the most general phase IVB. 
Hence we will sketch the derivation of the energy per 
particle for the latter. The elementary cell is a rhom¬ 
bus with the angle 4> between the vectors ai and 02 of 
the same magnitude a, see Fig. O The density of par¬ 
ticles on one plate is cr = I/(arsine/). We will prefer 
the parametrization of the angle by d = tan(0/2). An¬ 
other free parameter is a G [1/3,1/2] measuring the di¬ 
agonal shift c of the lattice on the opposite layer, see 
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FIG. 7. An excerption of the phase diagram for Yukawa par¬ 
ticles for the second-order phase transitions II —^ III and 
III —>■ IVA. The solid lines denote the critical lines obtained 
numerically by using Eq. ISSJ and (SH), respectively. The 
dash-dotted lines correspond to the asymptotic large-A for¬ 
mulas ISSli and (BZl). 


formula©. The square of the lattice vector can be writ¬ 
ten as = a^[(j-|-A:)^cos^(^/2)-|-(j —/c)^sin^(^/2)]. 

Next we distinguish the cases when j -|- fc is an even or 
odd integer and go to the summation over new indices m 
and n; details of this technicality and of the next steps 
can be found in Sec. Ill of paper [a. The main differ¬ 
ence is that we get (n -|- a)^ and {n — 1/2 -|- a)^ instead 
of and (n — 1/2)^ for the interlayer contribution. Ap¬ 
plying the Poisson formula m creates additional factors 
exp(27rzna) and exp[27rm(a — 1/2)]. Reducing the sum¬ 
mation over {— 00 , 00 } to { 1 , 00 } turns these factors to 
2cos(27rnQ;) and 2cos[27rn(Q; — 1/2)], respectively. The 
final formula for the energy per particle of phase IVB 
reads 




^3/2 ( 0, + fsj + 23/2 (^0, + f/5 


+4 X! + (-1)^+'"] ^3/2 (^0, - 7rzi/2 ^0, 


l + (-lv 

A 2 


00 r 


2ry^7r^ 


+ 2 E 


i=i 

00 


cos(27rja)z3/2 [tt r, /2, + j d ) + Z 3/2 [tt r, /2, + j /d 




2'K'^rf 


+2 E |cos 27r/ ^ 3/2 /2, + (-l)^Z 3/2 + J | 

E |cos(27rja) -h cos 27rj 0 (-1)''| 23/2 

i,h—\ 


00 r 


+ 2 E 


277 ^ 


1=1 

00 r r \2 


^3/2 ( ^ ) + ^3/2 ( J /^ 


2r}^ 


4 E ^3/2 + -27rzi/2 (^^>0 


+ E I ^3/2 

j,k— — oo ^ 

A" 


[ 21 ? 


2^^ /2 + tO’ + o) +k S 


-2'kzxii |^^,i?^/ 2 j + 4 E ^3/2 


i.fe=i 


E i 

2?7^ ’ S 


i,fc=i 
■ ^ 3/2 

r' 

2 


iyV2 + i(j + a - 1/2)^ + (fc - l/2fS 


fc-i) d 


+ ’^^(l + eE- 


r 


(40) 


VI. TRANSITION BETWEEN PHASES III AND The case d = 1 or ^ = 7r/2 is the fixed point of the 

IVA transformation S 1/S and corresponds to the critical 


One can verify that for the structure IVA with a = 
1/2 the energy (|40ll possesses the symmetry S —>■ 1/5. 
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point between phases III and IVA. In full analogy with 
the transition between phases II and III, we parametrize 
6 = exp(—e) so that the energy of phase IVA becomes 
an even function of e. The expansion of the energy (1401) 
around the critical point (5 = 1 in powers of small e takes 
the form 

AIVA(e“^??, A) = Aiii(I,? 7, A) + /i2(?7, A)e^ 

+h4(r], X)e^ + .... (41) 

The explicit formula for ^2 in terms of the generalized 
Misra functions is presented in Appendix D and /14 is 
also at our disposal. The critical line between phases III 
and IVA is once again given by vanishing of the prefactor 

h2(77^A=) = 0, (42) 

see Figs. |4]and[7l 


A. Critical behavior 


The expansion of the coefficients /12 and around the 
critical point ( 77 '^, A°) is analogous to the previous case 
of the second-order transition between phases II and III. 
The leading terms are /i 2 (??, A) ~ h 2 i{XX){'q — 77 ^) and 
114 ( 77 , A) ~ h 4 o(A‘^), where h 2 i(A‘^) < 0 and h 4 o(A'^) > 0 
for all Ac. Optimizing the energy Aiva with respect to e, 
the stationary solution e* = I — behaves as 


e = - 


fe2(?7,A) 
2 / 74 ( 77 , A) y 


/t2l(A^) 

2/74o(A‘^) 


1/2 




(43) 

with 77 —>■ ( 77 °)+. The order parameter e* has again the 
singular behavior of mean-field type with critical index 
/3mf = 1/2. We tested this results numerically in a plot 
analogous to Fig. [5]and got the slope (3 ~ 0.499. Without 
going into details, also other critical indices attain their 
mean-field values ( 1 ^ . 


B. Coulomb A —7 0 limit of the critical line 

The week screening (small A) case of the phase tran¬ 
sitions IIFIVA and IVA-V was studied by Monte Carlo 
methods in Ref. [ 2 ^. 

We reproduce 77 ^( 0 ) = 0.6214809246 [l3| in the 
Coulomb A —>■ 0 limit. The asymptotic shape of the crit¬ 
ical line for small A is again parabolic, see Appendix D: 

A^ a: C 34 [ 7 ?“ - 7 ?''( 0 )], C 34 ~ 149.7837254. (44) 

This asymptotic result is compared with the numerical 
calculation of the critical line directly from the relation 
(l42l) in Fig. [8l 


C. Hard-spheres A —00 limit of the critical line 


In the hard-spheres limit A —>■ 00 , the critical point 
for the III —?> IVA transition is (rjF 1/2 the 



FIG. 8. Transition III-IVA near the Coulomb A —^ 0 limit. 
The full line follows from the numerical treatment of the rela¬ 
tion 172 ( 77 , A) = 0. The dashed line corresponds to the asymp¬ 
totic formula dm. 


same as in the previous case of the II —^ III transition. 
Let us analyze the large-A limit of the critical relation 
FivF) = 0 - In the same way as for 772 , we get three 
leading terms from the seventh and ninth (last) sums of 
Eq. (IDip : 

16^^/" VV ’ T 4j VV ’ ^ 4 j 

“ 2 ^®/^ (v’ 2 ) •^>l-(45) 

The application of the asymptotic relations (IA9I) to this 
equation implies 



_^+ 

.877(772-^!)^/^ V Ay/^ 2 ^ 


\ 

A2(772 + i) j 


I 




= 0.(46) 


The root of the expression in the largest parentheses 
yields 


V 


1 

72 




In A ^ 



51n2 

2A 

+ 0 



(47) 


This asymptotic formula, taken for 77 ^, is plotted in Fig. [7] 
by the dash-dotted line. The comparison with the numer¬ 
ical results for the critical line (solid line) is very good. 
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VII. PHASE TRANSITIONS IVA - V, IVA - IVB 
AND IVB - V 

All phase transitions IVA - V, IVA - IVB and IVB - V 
are of first order due to a discontinuous change of both 
structure parameters S and a. For phase IVB one has 
to minimize numerically the energy (1401) with respect to 
two parameters S and a, which is tedious but feasible. 

We found that for A < 27.4436 phase IVA goes over di¬ 
rectly to phase V without entering the intermediate phase 
IVB. On the transition line, the parameter a jumps from 
1/2 to 1/3. In the Coulomb limit A —>■ 0 we get ??*(0) = 
0.732416, (5* = 0.69334 for phase IVA [l^ whereas for 
the rigid phase V = tan( 7 r/ 6 ) = I/-\/3 ~ 0.57735. The 
shape of the transition line is again parabolic, we can 
approximate it empirically by 

A^ ~ C 4 A 5 [??‘( 0 ) - r]*], C 4 A 5 ~ 805.3, (48) 

but now the parabola is reversed giving rise to the mul¬ 
tiple reentrant behavior, see Figs. |4]and|9l 



FIG. 9. A detailed view of the phase diagram around the 
tricritical point. The sector of phase IVA is connected by a 
narrow channel. 


The non-trivial 6 * for phase IVA increases to approx¬ 
imately 0.926 at A = 14 and then slightly decreases to 
0.853081 at A = 27.4436. In Ref. [IJ it was anticipated 
that there exist two disjunct regions of phase IVA in the 
phase diagram. Our more precise calculations indicate 
that there exists a narrow connecting channel merging 
these two regions into one, see Fig. |9l The maximum 
value of 5* for phase IVA is achieved when the chan¬ 
nel is the most narrow so that it does not decrease too 
much from the value i5 = 1 for phase III. Looking at 
Figs, m and |9] we can confirm the double reentrant sce¬ 
nario IVA-V-IVA-III-IVA-IVB [ 2 ^, restricted to a more 
precise interval 0.5275 <rf < 0.53643. 

For A > 27.4436, the phase IVB takes place and we 
have first-order transitions IVA-IVB and IVB-V, see Figs. 
[Hand [SI 


As concerns the transition line IVA - IVB, for A —?> 00 
it should asymptotically approach the value ( 77 *)^ —» 1/2 
so that phase IVA is absent in the hard-spheres limit . 
The value of <5* in phase IVA increases from 0.85308 at 
A = 27.4436 towards I for very large A. Concerning phase 
IVB, (5* increases from 0.763284 at A = 27.4436 to I for 
very large A and the other parameter a* from 0.41358 to 
0.5 along the same transition line. Thus, in the hard- 
spheres limit, the values 5* —^ 1 and a* —>■ 1/2 of phase 
III (see the top of Fig. |3]) will be attained as expected. 

Still in the hard-spheres limit A —>■ 00 , the transition 
line IVB-V should reach the point if ~ 0.877... [^ . 
Numerically, we got mere rf = 0.864133... even for A = 
500. The convergence is rather slow again, we have 5* = 
0.58102 and a* = 0.334428 for phase IVB at the same 
A = 500 value, gradually approaching the values 0.57735 
and 1/3 of phase V, respectively, with 0{l/\) corrections 
of both structure parameters. Now we want to derive the 
above hard-spheres result from our formalism. We recall 
that the energy ifivB is given by Eq. (jiOl) and Ey is 
its special case for 6 = l/-\/3 and a = 1/3. We apply 
the asymptotic formulas (IA 8 I) and (1X91) to the A —>■ 00 
limit of Eq. (piU) and neglect exponentially small terms. 
Five summands remain dominant; one from the sixth sum 
with j = 1, three from the eighths (last but one) sum, 
namely both terms with j = k = 0 and the second one 
with J = 0 , /c = 1 plus the j = k = 0 term from the ninth 
(last) sum: 


Eivb 



2 ^ 3/2 



+ Z3/2 


+2 Z3/2 


+4 Z 3/2 


A^ 77 ^ 

2772’ 2 
V 2772 ’ A5 




1 / 2 )^ 

5 


V 

4 


(49) 


Notice that two identical terms merged to the one on the 
third line. All these summands should be of the same 
order for very large A. Since the asymptotic relations 
(IA9I) imply that Zv{x,y) oc exp(—2yaly) for x —>■ 00 , 
and the hrst argument x = /{ 2 rf ) is common for the 
summands, the second arguments must coincide as well. 
Thus we have 



This equalities yield the expected asymptotic values of 
the structure parameters 5 = l/vX and a = 1/3. Simul¬ 
taneously, 


2 


0.877383..., 


A —>■ 00 . (51) 


This value can be rederived from purely geometric con¬ 
siderations, too. We have already mentioned that 
the particle density at one plate in phase V is cr = 
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l/[a^ sin( 7 r/ 3 )] = 2/(a^-\/3). For dense packed hard 
spheres of radius a, the perpendic ular distance of two 
layers of triangular lattices is d = y^ 2/3 a, see e.g. [^ . 
Inserting these values into rj = d^fa yields immediately 

(ED. 

We confirm that in the hard-spheres limit the transi¬ 
tion IVB-V will undergo no stepwise changes of struc¬ 
ture parameters and it will be of the second order, as 
expected. 

In Fig. [TUI present the transition values of the 
structure parameters 5* and a* for phases IVA and IVB 
at first-order transitions IVA - V (left, A S [0,27.4436]) 
and IVB - V (right, A G [27.4436, ooj). The left and 
right line fragments are separated by a gap, illustrating 
the step-wise change of structure parameters when going 
from phase IVA to IVB. We recall that the parameters 
of phase V are always fixed to 5* = l/y/Z and a* = 1/3. 
We found a tricritical point at = 0.772814 and A'^ = 
27.4436 where the three phases IVA, IVB and V coexist. 



FIG. 10. The transition parameters 5* and along the phase 
transition lines IVA-V (left) and IVB-V (right). Along two 
line fragments, A increases from 0 to oo. For the left fragment, 
the parameter a = 1/2 for phase IVA and a = 1/3 for phase 
V. There is a discontinuity in the parameters S and a between 
phases IVA and IVB at the tricritical point with A = 27.4436. 
The values of a on the right fragment correspond to phase 
IVB. In the hard-spheres limit A —>■ oo, the line ends up at 
the critical point ( 77 = 0.877383, 5^ = l/%/3). 


VIII. THE ENERGY PLOT 

We want to compare the values of the optimized en¬ 
ergy per particle for various values of A and 77 . We plot 
E{rf) for several fixed values of A in Fig. [TTJ These en¬ 
ergies vary by orders of magnitude, thus we have chosen 
semilogarithmic scale. 

First we consider two limiting cases. For 77 ^ 1 and 
A > 0, according to ([Tsp and (|22D the energy of the cor¬ 
responding phase II ln(i?n) ~ — 3 ^/^A /?7 and so In A di¬ 
verges if 77 —>■ 0. More interesting is the optimal energy 


of phase V, Ay, for 77 1. We were used to get the 

Coulomb limit as A —^ 0, but we can obtain this limit 
also for medium A and very large 77 , as the ratio A /77 —>■ 0 
again. For 77 ^ A, using the asymptotic formulas for 
the generalized Misra functions (Appendix A) we obtain 
from (|40D that 


F/y — TT — (1 -h e ^) + Cm ^ E 0 { 1 ), ( 52 ) 

where cm = —1.9605158... is the Madelung constant of 
the Coulomb potential for the hexagonal lattice; for an 
explicit representation of the Madelung constant in terms 
of z^(0, y) functions, see Eq. (24) with A = and 77 = 0 
of Ref. 17| . The leading term is the (minus) background 
energy (USD. 



FIG. 11. The dependence of the energy per Yukawa par¬ 
ticle E on the dimensionless distance rj for four values of 
A = 1,10, 20, 35, in semilogarithmic scale. 



FIG. 12. The derivative dEjdr] for A = 20. Full circle cor¬ 
responds to the second-order transition III-IVA, dashed line 
marks the discontinuity at the first-order transition IVA-V. 

We see in Fig. [TT]for few fixed values of A that the en¬ 
ergy is a monotonously increasing function of the dimen¬ 
sionless distance 77 . This means that the force between 
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the plates is always attractive. The non-analyticities 
at transition points are not clearly manifested in this 
scale. Therefore, for A = 20, we performed the deriva¬ 
tive dE/dr], directly for rigid structures and numerically 
using Ely A minimized with respect to 6 for phase IVA. 
The obtained results are plotted in Fig. [T2l We see 
the expected continuous but non-analytic behavior at the 
second-order transition point III-IVA as well as a jump 
discontinuity at the first order transition IVA-V. 


IX. INTERNAL PARAMETERS OF THE 
PHASES NEAR HARD-SPHERES LIMIT 

In and close to the limit of hard spheres A —>■ oo, the ex¬ 
pressions for the energies of the structures in terms of the 
generalized Misra functions admit an asymptotic analy¬ 
sis. This fact permits us to determine the ? 7 -dependence 
of the structure parameters of the present soft phases II 
and IVB in the A —>■ oo limit and eventually to derive 
their leading correction for large but finite A. 


A. Aspect ratio A of phase H at and near hard 
spheres 


The dependence of the aspect ratio Aps on rj for 
II is well known in the hard-spheres limit A —>■ oo 


phase 

m 


Ahs(? 7) = v '4774 + 3 - 2?7". (53) 


In the following, we derive this result and the first I/A 
correction to it by using our method. 

For A ^ I, most of terms in the energy of phase II 
m become exponentially small (we exclude from the 
discussion trivial terms which do not depend on A); only 
the term j = 1 in the sixth sum and the term j = k = 1 
in the eighth (last) sum contribute. As soon as A > 1, 
using (1X911 we get 


Eu 


V^A 


Z3 


A 2 1 


4?7^ ’ A 


—— 1 2z3 


1 

A 


V^Ae 


A^ 2 A 1 

+T+4A 

_A./r72_L A , 

T]V 'I ^ 4 ^ 4A 


+ f + 43 


ze ^ 


(54) 


The minimum of the energy is given by dEn/dlS. = 0, 
which implies 


1 A 
-k 


2VA 2Ary 
I 


0 77\/A = 


1 - ^ 




4772 -h A -f A 


+ f + IS ^ 


-j_ — I e“'^'\/''^+T'+43 


■ (55) 


significant than rationals and their arguments must be¬ 
come the same, i.e. \/rf E A/4 -|- 1/(4A) = 1 /\fK which 
leads to the known result dSSl). Numerics suggests that 
the correction is of the type 1/A, i. e. 


+ = + (56) 

We put the exponentials on one side, insert (l56l) and ex¬ 
pand yjif -I- A/4 -I- 1/(4A) — l/yfK up to the order 1/A. 
The absolute term vanishes and we have 


exp 


0 ( 77 ) 3 -I- 477 — 2rf' sj 477 "* -|- 3 
4?7 


„2 , A , J_ 
'/ + 4 - 1 - 4 ^ 


(\/4?7‘* -I- 3 — 2rfY/‘^ 
1 


_ /'57'\ 

(i “ 4 ^) 1-I -4774 - 2772 ^ 4774 ^ 3 ’ 

where we considered A Ri Ajjs on the second line. From 
this relation we readily get 

3/2 


7 ( 77 ) = 


(•\/4??4 -I- 3 - 2rf-^ 477 

3-^4774 - 2 Tfy/E^fTi 


X In 


1 -I- 477 "^ — 277^-^/4774 - 1 - 3 


(58) 


The value of 0 ( 77 ) is negative in the whole interval 0 < 
77 < \/y/2 of phase II. We find 0 ( 77 ) Ri —8 x 3 ^/^ 77 ^ for 
77 <C 1, confirming once more that phase II is entered 
directly from phase I for any small positive 77 . We tested 
the asymptotic result (l56ll . (l5^ numerically, see Fig. [TH 



FIG. 13. The aspect ratio A of phase II vs. A for four values 
of ?7 = 0.1, 0.3, 0.4, 0.5. The solid lines correspond to numer¬ 
ical calculations. The asymptotic X ^ 00 result (1561) . (1581) is 
represented by dashed lines. 


B. Parameters 5 and a of phase IVB in the 
hard-spheres limit 


If we want to reproduce just the hard-spheres limit As was shown above, in the hard-spheres limit A —^ 00 
A —>■ 00 , we can say that exponentials are by far more phase IVB takes place in the interval 77 € [l/-\/ 2 , 2 / 3 ^/ 4 ], 
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Let us study the A —>■ oo limit of the energy i?ivB dini); 
terms which do not depend on ^ or a are automatically 
excluded from the discussion. For very large A and gen¬ 
eral r] only the last two sums contribute, the remaining 
sums are exponentially small. The eighth sum has three 
important terms - one from the first Z3/2 summand with 
j = k = 0 and two identical ones from the second ^3/2 
summand with j = k = 0 and j = 0 , fc = 1 . From the 
ninth (last) sum we take the j = k = 0 term. The result 
is 


Ewb 


2v^A 


Z3/2 


A^ 2 /r, 


2 z. 


'3/2 




3/2 


2772’^ ' + 5 


4 


( 59 ) 


We notice that two more terms can become important 
in special limits. The first is the j = — 1 , k = 0 term 
from the eighth sum, the first Z3/2 summand, which con¬ 
tributes only in the limit ^l(i. e. 77^—^ 1 / 2 ). The 
other possibly important term can be found in Eq. (Hi) 
as the first one in the bracket, but it plays role only if 
7^2 —4/ ( 3 \/ 3 ) and it can be omitted for general 77 as well. 

Now we apply the asymptotic relations (IA 9 I) to the en¬ 
ergy (1^. The optimization of the energy with respect 
to parameters 6 and a leads to the equations dE/d6 = 0 
and dEjda = 0 . What we get are certain products of ra¬ 
tional functions and exponentials. To have a non-trivial 
solution in the A —>■ 00 limit, the dominant exponentials 
must have the same arguments which yields 



77 V 2 


(«-l/2)^ 

5 


J 

4 



This set set of equations can be readily rewritten as 


a = i (1 -I- (5^) , — a -I- = 0. (61) 

The quartic equation for S follows 

5"^ - 2(5^ -h Srf5 -3 = 0. (62) 

The discriminant of this equation —2^2(3 _ 147^4 _|_ 2 'Jrf^ 
is negative for any 77. Consequently, we get two complex 
roots and two real ones. It turns out that one of the real 
roots is negative and the only physical - real positive - 
root is given by Cardano formulas as follows 


5 = -S[-q) + -. + 4 .+ 


8772 


( 63 ) 


where 


5 ( 77 ) = 


U 1 

r , ^ 32 1 

3 

- Q(n)\ 


( 64 ) 


with 

1 /3 

g( 77 ) = 2^/3 ^'2777^ - 7 -f 3 v' 9 -42774-P8I778) . ( 65 ) 

The value of a follows straightforwardly from the first of 
Eqs. (I6T|) . 

It is easy to check that the above formulas give the 
correct lattice parameters at the endpoints of the phase 
IVB region, namely we have (a = 1 / 2 , ^ = 1 ) at 77 = 
l/v^ (phase III) and {a = 1 / 3 , <5 = l/\/ 5 ) at 77 = 2 / 3 ^/"^ 
(phase V). We did not find in the literature the above 
specification of the structural parameters of phase IVB in 
the hard-spheres limit. The numerical test of the results 
for phase IV parameters is depicted in Fig. [Ml For a 
given 77, the dependence of the parameters <5 (top set) 
and a (bottom set) on 1 /A, obtained by numerics, is 
represented by open symbols (connected by solid line), 
the asymptotic A —^ 00 result given by our Eqs. den is 
depicted by full symbol. It is seen that numerical data 
converge quickly to their asymptotic values. 



FIG. 14. The structure parameters 5 and a of phase IVB vs. 
A for three values of 77 = 0.75, 0.8, 0.83. Numerical data are 
represented by open symbols (connected by solid lines), the 
asymptotic A —^ 00 result given by Eqs. m is depicted by 
full symbol. 


X. CONCLUSION 

In this paper, we have studied the zero-temperature 
phase diagram of bilayer Wigner crystals of Yukawa par¬ 
ticles. To calculate the energy per particle of the phases, 
we used the recent method of lattice summations da 
extended to Yukawa potentials. The weak point of the 
method is that one has to know ahead the possible phases 
from numerical simulations. The strong point is that the 
truncation of the series of the generalized Misra func¬ 
tions provides extremely precise estimates of the energy. 
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e.g. the truncation at the 5th term provides the accuracy 
within 17 decimal digits. 

Another strong point of Misra functions is that they 
can be readily expanded around the critical point, pro¬ 
viding in this way closed-form expressions for the critical 
lines between phases II-III (l25|) and III-IVA (I42|) . Only 
few Misra functions contribute in the equations for the 
critical lines in the asymptotic Coulomb A —>■ 0 and hard- 
spheres A —>■ oo limits. The characteristic feature of the 
Coulomb limit is the parabolic shape of the critical lines, 
see Eq. (1351) with the corresponding plot in Fig. [Blfor the 
II-III phase transition and Eq. (IHl) with the correspond¬ 
ing plot in Fig. | 8 ] for the III-IVA phase transition. In 
the hard-spheres limit, the asymptotic formulas for the 
II-III phase transition (IMl) and the III-IVA phase tran¬ 
sition (H71) are pictures by dash-dotted lines in Fig. [T] 
It turns out that the second-order phase transitions II- 
III and III-IVA exhibit the mean-field critical exponents 

dMl). 

The most important features of the Yukawa phase di¬ 
agram obtained by Messina and Lowen were con¬ 
firmed. On the contrary to previous suggestions, phase 
I goes directly to phase II at ?; = 0, i.e. there does not 
exist a finite interval of positive ry-values where phase I 


dominates. Another important novelty is that instead 
of the suggested region of phase coexistence, we found 
a narrow channel within one continuous region of phase 
IVA. This fact also lead to the tricritical point where the 
phases IVA, IVB and V meet. 

Another application of our formalism is the determina¬ 
tion of the structure parameters of soft phases II and IVB 
in and close to the hard-spheres limit. For A —> oo, the t]- 
dependence of the aspect ratio A of phase II has already 
been known see Eq. (15^ . We were able to derive 
the first 1/A correction to this asymptotic relation, see 
Eqs. (1551) and (1551) . which is in perfect agreement with 
the numerical results (Fig. 1131) . The derivation of the t]- 
dependence of two structure parameters S and a of phase 
IVB in the limit A —>■ oo, see the relations (1501) and Fig. 
m is likely new as well. 

As concerns future perspectives to apply our method 
to other systems, the system of particles with 1 /r'^ inter¬ 
actions seems to be a good candidate. 
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Appendix A 

We give explicit analytic formulas for several z^{x,y) functions p7l) with half-integer arguments: 


z.dx.v) = 


1^5 “*(vf ^ 

1 - t ./wy"] + erfc(Y^+ 




y 


2^J 4y3/2 _ 




+e-2%/^(l + 2yi^) erfc(^y^-+e2^(-l + 2yi^) erfc(^y^- 


■ 


' 




-4e-"/---^(3 + 2^y)V^ 


(Al) 

(A2) 

(A3) 


_l_g (Axy + + 5) (Axy — + 5) erfc^-I- 


.(A4) 


I- 

Here, we introduced the complementary error function the obvious relation 


2 

erfc(z) = —= / exp (—<^) dt. 

Jz 


(A5) 


dz^{x,y) 

dy 


= -z^+i{x,y). 


(A6) 


The case = 1/2 can be found at the end of Ref. [^ . 
The expressions for larger v can be obtained by applying 


The Misra function case z^{0, y) should be under- 
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stood in the sense of the limit a; —>■ 0 , 

^1/2(0,?/) = - -Ky/y erfc(V^)] , 

v’’" 

23/2(0,1/) = ^ erfc(v^), 

25 / 2 ( 0 , 1 /) = [2e~^y./y +erfc(V^)] 

27 / 2 ( 0 , y) = 2e-^yy/y (3 + 2TTy) 

+3 eT^c{^/^Ty) 


(A7) 


27 / 2 ( 3 ;, y) = 


Q-vy-x/TT 

y 


1 + Oi- 
y 


■ (AS) 


For y finite and x 1, we have 

23/2(3:, y) = + ^ . 


+ 0 




lAQi 


We need also the asymptotic expansions of z^{x,y) when 
one of the arguments x or j/ is large. For x finite and 
?/ > 1 , we get 


21/2(3;, y) 

Z3i2{x,y) 


e 


e 


-ny-x/ 

yj^3/2 

-■ny-x! 


TT 


TT 




Q-'^y-x! 

25/2(3:,y) = lA- 


1 + 0 1 - 

,y 


1 + 0 


1 + 0 



We applied the large-argument expansion of the error 
function [ 2 ^, erfc( 7 :) « exp(—z^)/(-Az). 

For small arguments 5x and Sy, we shall need the fol¬ 
lowing expansion 

(x + Sx,y + 6y) = / _Q-io^+Sx)t^-iy+Sv)/t 

Jo 

'^e-^^e-y/\l-5x t){l-5y/t) 

Ki z„{x,y) - Sx z^-i{x,y) - Sy z^+i{x,y), (AlO) 
where we kept only terms linear in small variables. 



Appendix B 

The coefficients /i, 2 (y. A) in Eq. (fTSl) are given by 


/i(y, A) = 


V 

2 AFA 


2 E 

. i=i 


J 25/2 I 0, 


+4 ^ “ T ) ^5/2 1 0 

j,k=l A 

00 

1=1 


477^772 
A 2 


+ /V3 - 


3 A’ 4 ^ 2^2 + ^ 


J 


4772772 ^3 




3 25/2 TT ?7 , 


+3^vi\-i 


2 


A' 3 

+ 


+4 ^(-A(-l)Mfc2-:(^)z5/ 

j,k—l ^ 


3 ^ 


A 2 


+ 


A2E 


7=1 


+ 4 ^ (*,»- 

7,fc=l 


^5/2 I Vd I - — 25/2 

\ / X 2 72 

2 y \ ^ 3 


3E'^/A47/2’^ 


17/2 

k^Vs 


3 V 4772 ’y3;j 


4 E 

/,fe=i 


(*- 1 / 2 )"- 


(j- 1 / 2 )" 


25/2 


A" 2 , (y — 1 / 2 )" , T /r,^2 /^ 

IV’" + ^ +(fc-l/ 2 ) 4 / 3 _ 


(Bl) 
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/2(??, A) = 


2 ^\ 


2E 

^ i=i 


oo r -4 
3 


— ^ 7/2 I 0 


18 


A^ f 

+ 




Y^7/2 ' 0 


’ Air'^rf 




A" 


47r2772 yay 


+ 4 E 

j,fe=i L 


1 


o ( “ "T ) ^7/2 ( 0 > 


A^ 


47r2?72 


+ 2 ^ + fc2y3 - 


A2 


E 

18 ' 


+2E(-1) 

1=1 

(22 

E" ’ 

00 

+4 ^ (-ir(-i) 

j,k=l 

f (22 
Zb/2 71 17 , 


2„2 


A2 


Z7j2 TT 77 , 


+ 


3 


A^ 


47r2772 yjz 
2 


+ -77Z7I2 1 TT T] 


- 2^2 


3^Z5/2 ^ 

A2 

’ 47r2?7' 


^ + k^V 3 


2 +f^ 


3 


47r2?72 y /3 

1 




3\/3 

00 

2E 

1=1 


A^ 


3 


18 


47r2772 ^ 

.4 


k'^Vs 


A^ 72 


'vV’TlJ 


A" .2 /:t\ 

4772 ’^ 3^3 


+ 4 E 

l.fe=i '■ 


fc 2-4 


^7/2 


X 2 42 \ -2 

^ J +fc2^E-^ 

477 ^ ^ 


Zb/2 


J.fc=l 

(j-1/2)" 

3 V 3 


(fc -1/2)2 - 


Z7/2 


3\/3 


rA 2 , , (j- 1 / 2)2 


4772 


73 


E E 

W' y/i) . 

a 2 72 

4772 ^3 

(fc- 1/2)273 


rA2 ^2 ^ O' - y ^ (fc _ ^^2)273' 


4772 


73 


(B 2 ) 


We are interested in the small-77 behavior of the above functions. One of the arguments in the Zi,{x,y) functions 
becomes large, thus we can apply the asymptotic relations (IA8I) and (|A 9 p . Neglecting the exponentially small terms 
we find that only the seventh and the ninth (last) sums both in (IBII) and (IB 2 p contribute, namely the leading terms 
with / = 1 and j = k = 1 , respectively. Consequently, for a fixed A > 0 and 77 —i 0 (i.e. A/77 1 ), we have 

A 2 1 A 2 / A 2 , 


/i(i 7 , A) 
/ 2 ( 7 ?, A) 


2-y7rA 


77 

'1 

2i/7rA 

9 

A 



Zb 12 


3 

77/2 


47 ’ 73 

A2 1 


75/2 


47 ’ 73 E 


4772 

A2 


,17 


8 3®/'^77 


e 3T7 s 


1 

EZJ 
TsZJ 


1 A 

0 31/^77 


273 


6 


( 7 ^ + 7 


(B 3 ) 


14477(7 + 

where we repeatedly neglected subleading terms. Expanding also the second exponential exp[—A 7 1 + 73772/(3^/‘*?7)] 
in 77, we get (1^ . 


Appendix C 


The coefficient (72(77, A) in Eq. (IMl) takes the form 
^ A 2 


92(77, A) = 


77 


T^A \ ^ 

'' \_7 = 1 L 


J‘ 2 : 7/2 I 0 , 


4772772 


+ J ) ~ 9 2 : 5/2 


’ 47r2?7' 


+ 7 


+ E - k^fzr/2 ^0, + j2 -P fc2^ - (j2 -P 7)05/2 ^0, + f + /c2^ 


1=1 


J Z 7 I 2 tt ?7 , 


4772772 


+ J - J 2:5/2 77 77 


-2.„2 


477277 ^ 


+ /' 
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+ £ if - fc')"^7/2 4^ +/ + fc") - if + fc")^5/2 J^+f + 


3,k=l 

oo 

-E 

i=i ■- 


J^^7/2 


47^2 ’- 


J - J Z^/2 




-E 

j,fc=i '- 


^ ( ^,7^ + ^^ ) - if + ^^)-5/2 ( + k 


477 ^ 


+ E 1 - ik - 1/2)^] 27/2 

^ A2 


r \2 


477 


,,77^ + (j-l/2)^ + (fc-l/2)^ 


~ [O' ~ 1/2)^ + (A: — 1 / 2 )^] Z 572 


^^,^^ + o-i/2r + (fc-i/20 


(Cl) 


Our next step is to analyze the small A behavior of the critical line between phases II-III which is given by 
52 ( 77 , A) = 0; here, we write (? 7 , A) instead of ( 77 °, A'^) to simplify the notation. There are two small quantities: A^ 
and 77 — 770 , where we denote the Coulomb transition distance 77 ‘^( 0 ) = 779 ~ 0.262760268246823 .... Applying formula 
(lAlOl). we demonstrate one soecific examule of the exoansion of the generalized Misra functions in dcll), up to terms 
linear in small variables A^ and 77 — 779 : 


^ 7/2 


2 2 
TT 77 


’ 477 ^ 77 ' 


+r 


2 : 7/2 ( 


2 2-2 
TT r]f),J 


)- 


477 ^ 77 , 


2^9/2 ( 7 r^? 7 o>/) - 271^779(77 - 779)2:5/2 (77^0, j^) . 


(C2) 


Here, we used that rf' = [779 + (?7 — 50 )]^ ~ Vo + 2779(77 — 779). The absolute terms, like the leading one on the r.h.s. of 
(IC2I) . are canceled by the definition of the critical point at A = 0: 52 ( 770 ,0) = 0. Thus we are left with 




+ c^2^\v - m) = 0, 


(C3) 


where 

00 

^{23) = [fzg/2 {0,f) -fzy2 {0,f)] /i^T^^Vo) 

1=1 

00 

- E ^if ~ i^^f + k^) - if + k^)z7/2 {0,f + fc2)] /( 47 r 2 ? 7 ^) 

j,k=l 

00 

[f^S/2 -fz7/2 ( 77 ^ 779 ,/)] 7(477^77^) 

i=i 

00 

- E [if - k^fzQ/2 ( 77 ^ 50 ./ + kf - if + k^f7/2 (77^5o>J^ + kf] /i^^T^^Vo) 

j,k=l 

00 

“E [4^^572 (0>/) -fz3l2 {0,f)] /i^Vo) 

1=1 

00 

~ E [if ~ f)‘^Zg /2 (0, + k^) — if + f)z3/2 { 0 ,f + fc^)] 7(4^9) 

j,k=l 

- E \ [(7 - 172)" -ik- 172)"]' Z 3/2 [0, 5 o + (j - 172)" + ik- 1/2)2] 

j,k=l ^ 

+ [(j - 1/2)2 ^ _ ^/2)2] ^ 3/2 [0 ,772 + (j - 1/2)" + {k- 1/2)"] |/(477g) -0.04791591901052 (C4) 


[/2:5/2 (77"77o,j") -f^3l2 (77^5o:/4] ‘^f'^0 


i=i 


and 


















18 


- [(f - {Tr^vl,f + k^) - if + k^)z3/2 + k"^)] 27r^»?o 

j,k=l 

- E \ [(■?■ - 1/2)" -{k- 1/2)2] 2 ^2 ^ _ 1/2)2 ^ _ 1/2)2] 

j,k=l k 

+ [{j - 1/2)2 (fc_ 1/2)2] _ 1/2)2 1/2)2] ']^2r,o fv 1.15830861669576. (C5) 

Eq. (IC3I) with the specified constants yields (1351) . 


Appendix D 


The function /i 2 (??, A) in Eq. (HTl) reads as follows 


^ 2 ( 77 , A) = 


E[i + (-i)'] 


00 

+ E [i + (-i)^'(-i)l 

j,k=l 


( 0 , ^ - fz,/2 ( 0 , ^ + f 

(j2 - k'^ fzT/2 f0, —^ + /2 + kA - [f + A:2)z 5/2 f0, —+ f + 




2TT^ri‘^ 


+E[i + (-i)^] 


i=i 


j^zr/2 ( r; 2 ^ 2 / 2 ^ + ^- 2 ^ _ ^ 2 ^^^^ (^^^^^2, + f 


^ r /22\2 \ / 

+2 ^ (-1)-^ {f - k^fz'rj2 27 t^ J ~ ^^)^5/2 

00 

+E 


22 \ 2 

7 ] 7 T X 


+f + e 


2 ’ 27:"^ rf 


1=1 ■- 


4 / -2 \ -2 

J Z7/2 { ^,J ) - J ^5/2 


2772 ’•■ 


+ E 

3 ,k=l 


(j" - k‘f‘v2 (- if + k)--,,, + k 


+ E i [(j - 1/2)"-(fc-1/2)2]" Z7/2 

l,fc=i 1 

~ [(/ ~ 1 / 2 )" + {k — 1 / 2 ) 2 ] 


2r] 


(j - 1 / 2 ) 2 + (fc- 1 / 2 )^ 


,, (j - 1/2)2+ (fc-1/2)2 


+2 E 1 [(/ “ 1/2)^ “ 2:7/2 


2r] 

+O'- 1 / 2 )^ + ^^' 


“ [(/ “ 1 / 2 )" + ^"] Z 5/2 


„2 

ilZa ’ T + (/ ~ 1 / 2 )" + k 


2 , 7.2 


+ E 1 (/ “ 1 / 2 )"^ 2 ; 7/2 

i=i 


2772 ’ 2 
r \2 ^2 


^ + (j ~ 1 / 2 )" — (j — 1/2)2 z5/2 


X 2 T 72 


(Dl) 


Now we can analyze the low-A limit of the critical line between phases III and IVA. Proceeding in the same way as 
in the previous case of the II-III transition, taking the value 770 = 77 '^( 0 ) ~ 0.621480924579783, we get from (ID1|) the 
equality 


cf'^^A2 + cE(7?-77o) = 0 

with ft! 0.0063328359292865 and si -0.94855575801235884369, so that Eq. ^ follows. 


(D2) 
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